DD-CSA cross-correlation in a two-spin system

Ilya Kuprov,
University of Oxford.
kuprov@physchem.ox.ac.uk

This is a simple program for calculation of diamagnetic dipolar and CSA relaxation, DD cross-relaxation and DD-CSA cross-correlation rates between user-specified nuclear spin orders in a two-spin system. The internuclear distance is assumed to be fixed (rigid molecule). Both dipolar and chemical shielding tensors are assumed to be axial, an assumption that is strictly valid for inter-nuclear dipolar interaction, but may sometimes be violated for shielding tensors. Specify the source and the target spin order in the program body and evaluate the entire worksheet. β is  an angle between the internuclear vector and the axis of the shielding tensor. Pencil evaluation of these relaxation rates is very bulky, the best places to look are:

M. Goldman, Interference effects in the relaxation of a pair of unlike spin-1/2 nuclei, J.Magn. Res., 60 (1984), 437-452.
M. Goldman, Formal theory of spin-lattice relaxation, J. Magn. Res., 149 (2001), 160-187.
A.G. Palmer III, in: Protein NMR spectroscopy: principles and practice, Academic Press, 1996.
M. Mehring, V.A. Weberruss, Object-Oriented Magnetic Resonance, Academic Press, 2001.

The calculation uses group theoretical properties of the interaction Hamiltonian, expanding it into a sum of irreducible tensor operators, which are rotated using Wigner matrices. The Wigner matrix elements are then autocorrelated in the Bloch-Redfield-Wangsness theory treatment.

General setup:

Needs["LinearAlgebra`MatrixManipulation`"]; 
      Off[General::spell1, ClebschGordan::phy];
Id=(({{1, 0}, {0, 1}})); (Subscript[σ, x])=(({{0, 1/2}, {1/2, 0}})); (Subscript[σ, y])=(({{0, -I/2}, {I/2, 0}})); (Subscript[σ, z])=(({{1/2, 0}, {0, -1/2}}));
Kron[a_, b_] := BlockMatrix[Outer[Times, a, b]];Comm[A_, B_] := A . B - B . A;Scal[A_, B_] := Tr[A . Conjugate[Transpose[B]]];
(Subscript[H, x])=Kron[(Subscript[σ, x]), Id];(Subscript[H, y]) = Kron[(Subscript[σ, y]), Id];(Subscript[H, z]) = Kron[(Subscript[σ, z]), Id];(Subscript[H, plus]) = (Subscript[H, x]) + I (Subscript[H, y]);(Subscript[H, minus]) = (Subscript[H, x]) - I (Subscript[H, y]);
(Subscript[F, x]) = Kron[Id, (Subscript[σ, x])];(Subscript[F, y]) = Kron[Id, (Subscript[σ, y])];(Subscript[F, z]) = Kron[Id, (Subscript[σ, z])];(Subscript[F, plus])=(Subscript[F, x]) + I (Subscript[F, y]);(Subscript[F, minus]) = (Subscript[F, x]) - I (Subscript[F, y]);
(Subscript[DQ, x]) = 
      (1/(Sqrt[2])) (2 (Subscript[H, x]) . (Subscript[F, x]) - 2 (Subscript[H, y]) . (Subscript[F, y])); 
    (Subscript[DQ, y]) = (1/(Sqrt[2])) (2 (Subscript[H, y]) . (Subscript[F, x]) + 2 (Subscript[H, x]) . (Subscript[F, y])); 
    (Subscript[ZQ, x]) = (1/(Sqrt[2])) (2 (Subscript[H, x]) . (Subscript[F, x]) + 2 (Subscript[H, y]) . (Subscript[F, y])); 
    (Subscript[ZQ, y]) = (1/(Sqrt[2])) (2 (Subscript[H, y]) . (Subscript[F, x]) - 
            2 (Subscript[H, x]) . (Subscript[F, y]));
WignerD=(({{Cos[(β/2)]^4, (1/2) Sin[β] (1 + 
                      Cos[β]), (Sqrt[(3/8)]) Sin[β]^2, -(1/2) Sin[β] (Cos[β] - 
                      1), Sin[(β/2)]^4}, {-(1/2) Sin[β] (1 + 
                      Cos[β]), (1/2) (2 Cos[β] - 
                      1) (Cos[β] + 
                      1), (Sqrt[(3/2)]) Sin[β] Cos[β], (1/2) (2 Cos[β] + 1) (1 - 
                      Cos[β]), -(1/2) Sin[β] (Cos[β] - 1)}, {(Sqrt[(3/8)]) Sin[β]^2, -(Sqrt[(3/2)]) Sin[β] Cos[β], (1/2) (3 Cos[β]^2 - 
                      1), (Sqrt[(3/2)]) Sin[β] Cos[β], (Sqrt[(3/8)]) Sin[β]^2}, {(1/2) Sin[β] (Cos[β] - 
                      1), (1/2) (2 Cos[β] + 1) (1 - 
                      Cos[β]), -(Sqrt[(3/2)]) Sin[β] Cos[β], (1/2) (2 Cos[β] - 1) (Cos[β] + 
                      1), (1/2) Sin[β] (1 + Cos[β])}, {Sin[(β/2)]^4, (1/2) Sin[β] (Cos[β] - 
                      1), (Sqrt[(3/8)]) Sin[β]^2, -(1/2) Sin[β] (1 + Cos[β]), Cos[(β/2)]^4}}));
$Assumptions = (Subscript[ω, 1]) ∈ 
            Reals && (Subscript[ω, 2]) ∈ Reals && τc > 0 && 
          Dip ∈ Reals && ℏ > 0 && 
          CSA ∈ Reals && β ∈ Reals;

(Subscript[T, 2, -2]) = 
        (1/2) (Subscript[H, minus]) . (Subscript[F, minus]);
(Subscript[T, 2, -1]) = 
        (1/2) ((Subscript[H, z]) . (Subscript[F, minus]) + 
              (Subscript[H, minus]) . (Subscript[F, z]));
(Subscript[T, 2, 0]) = (Sqrt[(2/3)]) ((Subscript[H, z]) . (Subscript[F, z])  - 
              (1/4) ((Subscript[H, plus]) . (Subscript[F, minus]) + 
                    (Subscript[H, minus]) . (Subscript[F, plus])));
(Subscript[T, 2, 1]) = -(1/2) ((Subscript[H, z]) . (Subscript[F, plus]) + 
              (Subscript[H, plus]) . (Subscript[F, z]));
(Subscript[T, 2, 2]) = (1/2) (Subscript[H, plus]) . (Subscript[F, plus]);

(Subscript[A, 2, -2]) = 0;
(Subscript[A, 2, -1]) = (1/2) (Subscript[F, minus]);
(Subscript[A, 2, 0]) = (Sqrt[(2/3)]) (Subscript[F, z]);
(Subscript[A, 2, 1]) = -(1/2) (Subscript[F, plus]);
(Subscript[A, 2, 2]) = 0;

Shorthands:

\(ω\_1 = \(-γ\_1\)\ B\_0, \ \ \ \ \ \ \ \ \ \ \ \ ω\_2 = \(-γ\_2\)\ B\_0, \ \ \ \ \ \ \ \ \ Dip = \(γ\_1\ γ\_2\ ℏ\ μ\_0\)\/\(4\ π\ r\^3\), \ \ \ \ \ \ \ CSA = \(γ\_2\ B\_0\ Δσ\)\/3\)

\(Source\ and\ target\ spin\ orders\ in\ the\ laboratory\ \(\(frame\)\(:\)\)\)

Source = (Subscript[F, z]);
Target = 2  (Subscript[H, z]) . (Subscript[F, z]);

Bloch - Redfield - 
          Wangsness theory processing:

Hst = (Subscript[ω, 1]) (Subscript[H, z]) + (Subscript[ω, 2]) (Subscript[F, z]) // Simplify;
Hdip=Total[((({{(Subscript[D, -2, -2])[
                                  t], (Subscript[D, -2, -1])[
                                  t], (Subscript[D, -2, 0])[
                                  t], (Subscript[D, -2, 1])[
                                  t], (Subscript[D, -2, 2])[t]}, {(Subscript[D, -1, -2])[
                                  t], (Subscript[D, -1, -1])[
                                  t], (Subscript[D, -1, 0])[
                                  t], (Subscript[D, -1, 1])[
                                  t], (Subscript[D, -1, 2])[t]}, {(Subscript[D, 0, -2])[t], (Subscript[D, 0, -1])[
                                  t], (Subscript[D, 0, 0])[t], (Subscript[D, 0, 1])[
                                  t], (Subscript[D, 0, 2])[t]}, {(Subscript[D, 1, -2])[t], (Subscript[D, 1, -1])[
                                  t], (Subscript[D, 1, 0])[t], (Subscript[D, 1, 1])[
                                  t], (Subscript[D, 1, 2])[t]}, {(Subscript[D, 2, -2])[t], (Subscript[D, 2, -1])[
                                  t], (Subscript[D, 2, 0])[t], (Subscript[D, 2, 1])[
                                  t], (Subscript[D, 2, 2])[t]}})).(({{0}, {0}, {-Dip (Sqrt[6])}, {0}, {0}})))*(({{(Subscript[T, 2, -2])}, {(Subscript[T, 2, -1])}, {(Subscript[T, 2, 0])}, {(Subscript[T, 2, 1])}, {(Subscript[T, 2, 2])}})),2]//Simplify;
Hcsa=Total[((({{(Subscript[D, -2, -2])[
                                t], (Subscript[D, -2, -1])[
                                t], (Subscript[D, -2, 0])[
                                t], (Subscript[D, -2, 1])[
                                t], (Subscript[D, -2, 2])[t]}, {(Subscript[D, -1, -2])[
                                t], (Subscript[D, -1, -1])[
                                t], (Subscript[D, -1, 0])[
                                t], (Subscript[D, -1, 1])[
                                t], (Subscript[D, -1, 2])[t]}, {(Subscript[D, 0, -2])[t], (Subscript[D, 0, -1])[
                                t], (Subscript[D, 0, 0])[t], (Subscript[D, 0, 1])[
                                t], (Subscript[D, 0, 2])[t]}, {(Subscript[D, 1, -2])[t], (Subscript[D, 1, -1])[
                                t], (Subscript[D, 1, 0])[t], (Subscript[D, 1, 1])[
                                t], (Subscript[D, 1, 2])[t]}, {(Subscript[D, 2, -2])[t], (Subscript[D, 2, -1])[
                                t], (Subscript[D, 2, 0])[t], (Subscript[D, 2, 1])[
                                t], (Subscript[D, 2, 2])[t]}})).(WignerD.(({{0}, {0}, {-CSA (Sqrt[6])}, {0}, {0}}))))*(({{(Subscript[A, 2, -2])}, {(Subscript[A, 2, -1])}, {(Subscript[A, 2, 0])}, {(Subscript[A, 2, 1])}, {(Subscript[A, 2, 2])}})),2];
Hdn[t_] = 
            Hdip + Simplify[Hcsa, Trig -> False];
Dcomm = -Comm[Hdn[t], 
                Comm[Conjugate[
                    Transpose[
                      MatrixExp[I Hst τ] . 
                        Hdn[t + τ] . 
                        MatrixExp[-I Hst τ]]], 
                  Source]];
Rate = 
          Scal[Dcomm, Target]/(Sqrt[Scal[Target, Target] Scal[Source, Source]]); Rate = Simplify[Rate, Trig -> False] // ExpandAll;
For[a = -2, a ≤ 2, 
          For[b = -2, b ≤ 2, 
            For[c = -2, c ≤ 2, 
              For[d = -2, 
                d ≤ 2, 
Rate = 
                  ReplaceRepeated[
                    Rate, {(Subscript[D, a, b])[t] Conjugate[
                            (Subscript[D, c, d])[t + τ]] -> 
                        KroneckerDelta[a, c] KroneckerDelta[b, 
                            d] (1/5) (E^(-τ/τc))}]; 
d++;]; c++;]; b++;]; a++;]; Rate = Simplify[Rate, Trig -> False];
Clear[a, b, c, d];
Rate = 
            Integrate[Rate, {τ, 0, ∞}, 
              GenerateConditions -> 
                False] ; Rate = ComplexExpand[Rate] // Re // 
              Simplify ;

Source to target cross - relaxation rate:

Rate